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Expressions are developed for the transmission and reflection coefficients for propagation 
of a plane wave through a layered medium, taking account of the effects of the static 
magnetic field. A matrix formulation is used which allows proceeding to the limit of a 
continuously varying medium, and series expansions of the fields for this case are developed. 
The results arc expected to have application to interpretation of VLF data obtained within 
and above the lower ionosphere. 

1. Introduction 

Theoretical investigations of VLF propagation have primarily been concerned with the 
region below the ionosphere, as evidenced by the use of semi-infinite model ionospheres not 
only in earlier semiquantitative studies [Budden, 1951] but also as a reasonable simplifying 
assumption in current studies [Wait and Walters, 1963a, 1963b]. Even in studies in which 
a semi-infinite ionosphere is not assumed, the transmission coefficients are often not developed 
[Wait, 1962a, 1963] or are developed as an incidental part of the analysis [Johler and Harper, 
1962]. This emphasis lias been consistent with the nature of the available experimental 
data. However, as data from the ionosphere become available [Rorden et al., 1962], their 
interpretation requires consideration of propagation into and through the ionosphere. It is 
the purpose of this paper to provide a framework for this type of calculation. 

As is generally the case if tractability is to be maintained, some simplification of the 
physical situation is desirable. A suitable model for our purposes is a stratified plane ionosphere 
with plane waves obliquely incident. The earth's magnetic field is included since it is expected 
to allow propagation through the ionosphere (whistler modes) at some frequencies and angles 
of incidence. It is assumed to have constant strength but to be at an arbitrary angle to the 
plane of stratification. 

The equations are first developed for an arbitrary number of homogeneous layers, a 
problem considered by Johler and Harper [1962]. A different system of representation is 
used in this study, allowing the number of layers to be increased by multiplying together 
more and more matrices of constant (4X4) size rather than by increasing the size of a single 
matrix. This form has the advantage of allowing use of available computer subroutines 
to handle the matrices and the use of slow access storage for matrices not currently being- 
operated upon. The values of the matrix elements are determined by the solution of the 
Booker quartic for each layer. The terminology used is that employed by Yabroff [1957] 
in treating a single planar boundary. 

In the next step, it is assumed that the layers become infinitely thin while becoming 
infinite in number, to represent a continuously varying medium. A series expansion for the 
transmission and reflection coefficients is then obtained. The reflection coefficients series 
have the form found by Wait [1962a] and by Heading [1963] by an iterative procedure. A 
physical interpretation of the series terms as representing multiple reflections within the 
varying medium is obtained by examining the development of the equations for the limiting 
conditions. 
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2. Layered Medium 

2.1. Equations Holding in Each Homogeneous Slab 

Waves propagating in a homogeneous plasma medium are governed by Maxwell's equa- 
tions and (neglecting the ions) by the equations of motion for the electrons: 



VXE= 



"Mo 






VXH = -7VeV- 



«o 



dE 



m 



dy 



sE— mv\— /*oe(VXH ; 
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(3) 



where E and H are the wave fields, V the electron velocity, N the electron density, e the electron 
charge, m the electron mass, v the classical electron collision frequency, e the dielectric constant 
of free space, and /x the permittivity of free space. If the energy dependence of the collision 
frequency, v, is considered, v is replaced by a complex, frequency dependent quantity, as 
discussed by Wait [1962b] and Johler and Harper [1962]. This substitution does not affect the 
the equations to be derived below. 

These equations assume that ion motion can be neglected and that the plasma is "temperate" 
[Allis et al., 1963], i.e., that the electron motion occurs over only a small fraction of a wave- 
length and that the electron thermal speed is much larger than the induced speed. These 
assumptions imply that the wave magnetic field H is much less than the earth's field H and 
hence can be dropped from the last term on the right in (3). 

The notation used follows Yabroff [1957]; a right-hand cartesian coordinate system is 
defined with x normal to the plane of stratification, so that N and v are functions of x only. 
The direction of a vector in this coordinate system is specified by the angles a and /3. Beta is 
the angle between the vector and a unit vector in the x direction; a is the angle between the 
projection of the vector in the y-z plane and a unit vector in the z direction, as shown in figure 1. 
H is taken to lie in the x-z plane at an angle f3 H to the x axis. The incident wave in free space 
is taken to be propagating with its wave normal given by angles a T and p r . 

The incident wave will be of the form 



F I =(F I0 ) exp i- (ct — zcos Q I — ?/sin /3 7 sin a, — 2 sin &cos a,) 



(4) 




Figure 1. Coordinate 
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where (F I0 ) is the complex magnitude of the field at the origin of the coordinate system (taken 
to be in free space). The reflected wave will be similar, with a R = ai] (3r=(it— (3 f ). In each 
homogeneous layer, for such an incidenl wave, the" waves will have the form 



Fi=(F iQ ) exj) i— (ct — xD—ys'm /3 7 sin ctj — z sin ft cos a r ) 



(5) 



where the continuity of the fields at the boundaries between homogeneous layers has been 
employed to specify the y and z variation. For such waves in a layer, (1), (2), and (3) ran be 
manipulated to eliminate H and V, giving a set of homogeneous equations in the components 
of E: 



'mix — a 2 mi2+a T D rriiz-\-a L D 

— mi 2 -\-a T D m 2 2—T) 2 —a 2 L m^-\-a L a T 
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(6) 



The particular form of (he matrix multiplying E in (6) is obtained by deriving the conductivity 
tensor relating V to E from (3), eliminating H from (1) and (2), and replacing NeV in the result- 
ing equation by the conductivity tensor multiplying E. A detailed derivation of the equivalent 
(for a different choice of coordinates) of (6) is given by Budden [1961]. 
The parameters m u and (/ are 
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in terms of 



and 



a=sin p f , (a) 
a L =a cos a T , (b) 
a T =a sin a T , (c) 

Ullf) 
fr = 7? (a) 

Jo 

h L =h cos fi H (b) 
h T =h sin f$ H (c). 



(8) 
(9) 



(10) 



Finally, the quantities /, / , and f H are the familiar wave frequency, plasma frequency, 
and gyrofrequency, respectively: 



h- 



47r 2 e m 
2tviu 



(12) 



Nontrivial solutions to (6) exist only if the determinant of the matrix multiplying* E 
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vanishes. This condition produces a quartic in D to be solved for each layer: 

a±D* + a z D s + a 2 D 2 +aiZ?+a =0 (13) 

where 

a 4 = mil (&) 

a 3 =2a L m 13 (b) 

^2 = m 33 a£ + m 22 a r + m, i (a 2 — m 33 — m 22 ) + mf 2 + m? 3 (c) 

a l = 2m n {al-\-a L a 2 T ) — 2a L m 22 m n (d) 

0o= N 33 — a 2 T ) (m 22 — a£) (mn— a- 2 ) + (m 33 — 2 r )ra? 2 — (m /2 — a£)m? 3 

+ 2m 2Z m n m l2 —(ala 2 T —m 2 2Z ) (m n — a 2 ) (e) (14) 

Thus there are generally four modes possible in the medium, two upgoing and two downgoing. 

Equations (13) and (14) are the Booker quartic [Booker, 1939] in the form used by Yabroff. 
He needed only two roots for D corresponding to upgoing waves since only a single boundary 
with a semi-infinite medium on either side was treated. For a number of boundaries, all four 
modes have to be included except in the final semi- infinite medium. 

The properties of the wave in the medium are contained in D. The wave has its phase 
normal in the direction 

a M =ai (a) 

^^teD^ (15) 

which is Sn ell's law. The phase velocity is 

F * = [(Re£) 2 +sin 2 /y 1/2 ' (16) 

The wave is seen to be an inhomogeneous plane wave since it is attenuated exponentially in 
the x direction at a rate of 

f Im (D) (17) 

nepers per unit distance. 

2.2. At the Boundaries 

At each boundary the tangential components of E and curl E have to be continuous. 
Thus the equations can be written 

E y -=E y+ (a) 

(VXE),_ = (VXE), + (b) 

(VXE) 2 _ = (VXE) 2+ (d) (18) 

where ( — ) and (+) refer to the field just below and just above the boundary, respectively. 
The equations may conveniently be written in terms of matrices 

A i E i .=A w E i+1+ (19) 

where E- and E+ are column vectors of independent (one for each mode) components of E in 
the two mediums i and i+1 forming the boundary, and A t and A i+ i are matrices, characteristic 
of each medium, relating the elements of E- and E+ 9 respectively, to E vy (VXE)^^, and 
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(VXE) 2 . Thus (a^nE-i is the y component of E- u {a i+ i) u E +2 is the y component of E+ 2 , 
(a t ) 2 iE-i is the y component of (VXE)_!, and so on. 

For the incident and reflected waves below the first boundary, the independent components 
of E- are conveniently chosen to be the components of E in (E\\) and perpendicular to (E±) the 
plane of propagation. Thus 

'(E ± )^ 



#-.= 



(#ll)i 

(Ex); 



(20) 



where the subscripts 1 and 2 refer to the incident and reflected waves, respectively. The 
matrix A a is, then, 



A 



cos olj sin a 7 cos /3 7 cos a r sin a 7 cos j8/ 

—sin « 7 cos /3 7 cos a 7 (cos 2 ft — sin 2 ft) sin a 7 cos T ft —cos a 7 

— sin c* 7 cos a 7 cos ft —sin a 7 cos a 7 cos ft 

^— cos aj cos ft sin a T (sin 2 ft— cos 2 ft) cos a 7 cos*ft sin a 7 



(21) 



In the plasma mediums a more convenient choice of independent variables is E y for each 
mode: 



The matrices A: are then 



A = 



«(, 



V. a «41 



l 
a, 2 . 



(Ey), 
(Ey), 

<(E V )U 



(22) 



1 



1 ^ 



(23) 



where 



(a,i) 2 j=— cos a 7 sin PjS^+DjRj (a) 
(a«)ay=fij (b) 

(^i)4j = — ^Dj+sin a 7 sin ft-S j (c) . 



(24) 



The quantities i?j and S? are the cof actor ratios of the matrix in (6), which relate E z and 
E x , respectively, to E y for the^th mode in the ith layer: 



B}= 



( — m 1 2 + a jJ?}) ( — m 23 +a L a T ) — [m 22 — (ffj) 2 — ai](m 13 +az#*) 



(25) 



(^23+^L«r)(mi3+^A-) — (— mi2+a r 2?j)[m 3 3— (#j) 2 — 4] 

Si __^ [m 22 — (v9j) 2 — a£][m 33 — (A-) 2 ~ a 2 r]— (^s^-a^X— m 23 +a L a T ) ( 

(m 2 3+aL«7')(m 1 3+a/>D}) — (— m ]2 +a 7 Z>j)[m 3 3— (D}) 2 — a%] 

It must be remembered that the m is also characteristic of the ith medium. 
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Finally, at the last boundary, the reflected modes are not present in the upper medium, 
so that 



E N4 .= 



(E v ) 2 





(27a) 



where the subscript N is the number of boundaries. The form of A N need not be changed, 
however. 

If the final semi-infinite medium is free space, a better choice for the transmitted modes 
would be 



E, 



(£|i). 
o 





(27b) 



with A N given by (21). 

To relate the incident and reflected waves at each boundary of a layer, another matrix, Aj 
is defined: 

E i .=L i E i+ 



where A* is diagonal 



(28) 



-i- di Z»i 
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1 











e c 














e c 3 






(29) 



L e c ' * 

and where d t is the thickness of the ith. layer. 

With the above notation, the incident and reflected waves at the first boundary can be 
related to the transmitted wave at the last boundary by a series of matrix multiplications. 



En+ — \An A N -\&n-\A n -iA n - 2 &n-2A n -2 . . . Ai&iAi A )£L - 
We define the product matrix, M, so that (30) is more concisely written: 

E N+ =ME - 

M^{A- N 1 A N .^ N _ l A^l l . . . A^Ar'Ao). 



where 



(30) 

(31) 

(32) 



Note that if a layer, i=a, is allowed to vanish by letting d a become zero, A a becomes the 
identity matrix and the product A a A a Aa l is also the identity matrix. 
By splitting £" _ into the incident and reflected modes 



E _=Er+E R 



(33) 
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where E r and E R have the form 





'(^±>r 




E,= 


mi 




(a) 


Er = 




(E-L)l 

v.CE||)2> 


(b) 


(31) can be rewritten to give solutions E N+ (=E T ) and E R in terms of the source E,\ 


M~ 


'^-^ 


*=#/. 



(34) 



By defining a matrix, E„ that combines EV and E R , (35) can he written 

NE S =E, 
where 2Vis seen to be of the form 



(35) 



(36) 



N-- 



M- 






<r 








1 


-1 


1 


-1 



(37) 



where the notation means that the left-hand half is identical with the left-hand of M l and the 
right-hand half is as shown; while E s is the column vector 



i.e. 



E S =E T +E B 

(E y ) 2 
{E L ) 2 



E 8 



(a) 



(b) 



(38) 



Having solved the quartic, (13) for the roots, D, in each medium, the elements of the 
matrices forming M, and hence N, can be evaluated. The transmitted and reflected waves 
then are given as solutions of the four linear equations (36). The ratios of (E ± ) 1 to (E±) 2 
and (£11)2, and of (i?||)i to CEj|) 2 and (E±) 2 give the reflection coefficients ±R±, ±R\\, \\R\\ and 
\\Rj_. When the second semi-infinite medium is free space the transmission coefficients can 
be similarly defined. 



ri6-79T— 04- 
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Once the transmitted and reflected waves have been determined the fields within the 
medium can be found. At some boundary within the medium, we have, in the ith layer 
just above the boundary, 



■tLi+ — A% (Ai-i^i-iAi-i 



A^At^AA- 



(39) 



where the reflected wave components of E - are now known, so that they need only be inserted. 
There is an interesting" relationship between (36) and the expressions developed by Johler 
and Harper [1962] and Wait [1963]. Equation (36) is equivalent to the matrix equation (33) 
of Jollier and Harper, but it represents the effect of many slabs by a product 4x4 matrices 
rather than by a single matrix which grows in rank. Wait uses a method of development for 
a stratified ionosphere similar to that of Johler and Harper in treating the wave propagation 
problem in spherical coordinates where the earth and ionosphere are considered the boundaries 
of a waveguide. In this case, the resulting equation does not involve an arbitrary incident 
plane wave. E, and E T are solved for simultaneously (the earth provides another boundary, 
which would give two more equations and hence 6x6 matrices in the formulation employed 
here), so that the determinant of the equivalent of the matrix (N—I) for the spherical coordi- 
nates would yield the waveguide model equation. 

2.3. Limit of a Continuously Varying Medium 

The development in section 2 is particularly suited to proceeding to the limit of a con- 
tinuously varying medium. Swift [1962] employed it to investigate the characteristics of the 
coupling coefficients appearing in the differential equations. In this study we are more 
interested in developing forms for the reflection and transmission coefficient. Consider the 
form of the matrix M (32) . It is seen to be a product of elements 

T^AAiAT 1 (40) 

for each layer, except the semi-infinite layers on each end. The terms of I\ have the form 



( T i) ik — (A •) j\ (A i) \ m (Ai ) mk . 



By noting that A* is diagonal, 



-i-dD) 



(Ai) jk =8 jk e c 



) jk u jk 

one summation of (41) can be eliminated, leaving 

(rOfcH^OyiGAr 1 )!** 



-i-dD) 



(41) 



(42) 



(43) 



where the d { are all taken the same. 

As the thickness of each layer is decreased, V t approaches the identity matrix. By ex- 
panding A, in a power series for each term, we get 



A,= 



(i^dDi) 



1-i-dD* 
c 



-i-dD* 3 + 

c 



l-i-dDi+ ... 



. (44) 
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I\ can be written 

r _/_ i | r / r ; + ( 7 ; 2 ,/J' 7 »_ . . . . (45) 

Note that the i-d factors arc taken to multiply each term of the y t . The 7, have the form 

{y' i ) jk =(A i ) jl {AT l )iM (a) 

(t;%=(^)*(^%(A*) 2 (b) (46) 

and so on. 

The product of the matrices Y L can then be written 

n r,=[/-i|dyi+(i|d)V-. • .][/-i^dyi+. ..]... [i-i^+. . .]• (47) 

This can now be expanded: 

i = l C i = l \ C / i = l ; = i \ O / ? = 1 \ C / i = l jf=j 

i 7 '/) 2 2 tftf+( i ^ </) SSS tWK' + ... + ... (48) 

In the limit as d approaches zero, only the terms in 7' of (48) remain finite, there being on 
the order of (l/d) n terms in the summations for the nth term. For the terms containing 7" or 
a higher order, the number of terms in the summations increases as a lower power of (l/d) 
than the power of d appearing as a multiplier. Thus we are left with the series 

N OJ N / 0) V N N 

nr^/-i;r/s 7 ;+i^/ SS7W+... (49) 

i = l C j=i \ C J i=\ j=i 

as (/ becomes small. In the limit the sums become integrals, and we have 

M=A?-[l-i^f\'(x)dx+(i fffcdxrr'ix) f* dyy'(y)+. . ^A L (50) 

where A„ and A L are the matrices for the upper and lower semi-infinite mediums, respectively. 
This form of M can be evaluated If the quartic for D is solved with x a parameter to give a 
D(x) for each mode. 

The fields within the varying medium are now given by the limit of (39), which is 

E{x)=A-\x) [/-(/ £)£* y'(y)dy+(i ~J j^ dyy'(y)f* dzy'(z)-. . .] A L E (51) 

where the reflected wave elements of E Q are now known. 

The form for M given in (50) should be compared with the solution for the reflection 
coefficients given by Wait [1962a] and by Heading [1963] for a continuously varying medium. 
It is seen that the terms for the reflection and transmission coefficients in the present series 
have the same form as those produced by the iterative procedure used by Wait to solve a 
differential equation for the reflection coefficient and by Heading to solve an integral equation 
for the fields for an isotropic ionosphere. The present derivation is thought to give additional 
insight into the nature of the expansion, as well as to provide expressions for the transmission 
coefficients for an anisotropic ionosphere. Heading also derives an integral equation for the 
fields in an anisotropic ionosphere but does not develop the series solution. 

From the series in (48) it is seen that the higher terms express multiple reflections within 
the varying medium. This interpretation can be made clearer by considering a single boundary, 

415 



for which the equation for finite slab thicknesses was written 

A i E i .=A i+1 E i+1+ . (19) 

2?i+i_, the transmitted incident wave at the next boundary, is given by inserting (28) into (19) 
and rearranging: 

E i+1 -=(A i+l A^ 1 A i )E i _. (52) 

It is seen that E i+ i_ is given by the product of a matrix 1 of the form of T { with E t -. 

Consider now the relationship between the incident, reflected, and transmitted waves for 
an inhomogeneous region: 

E N+ =ME - (32) 

where E N+ is the transmitted wave, and E _ is the incident and reflected waves. Keeping (52) 
in mind, we can interpret the form of M given by (50) as follows. Using only the first term of 
the series in (50), we get the same relation as for a sharp boundary between the upper and lower 
semi-infinite mediums: 

E%+=A^A L E,_. (53) 

The second term alters this result by the reflection of the undiminished incident wave 
throughout the medium: 

E N+ = \_A-' (-; - c ^ U y'{x)dx} A L ~J E _. (54) 

The Az 1 and A Ly of course, couple this reflected wave to the semi-infinite mediums on either 
side of the inhomogeneous region. 

The third term then reflects the waA^e given by (54) throughout the medium 

E% + =A?{ -i - c £" dxy'(x) [£*" dyy'iyiJlArMo- (55) 

where now the inner integral from x to x u gives the wave at x that is due to the incident wave, 
E _, being reflected from the region beyond x; the outer integral then gives the result of reflec- 
tion of this wave throughout the medium. Higher terms in the series are similarly interpreted. 
This interpretation of the series terms agrees with that given by Westcott [1962], who considered 
the scattering process from infinitesimal slabs of the medium. 

We can check this for a single homogeneous slab, in which D(x) and, hence, y' are constant. 
Doing the indicated integrals in (48) produces the series 

M=A- u i{l-i^y'(x u -x L )-(~Jy' 2 ^^)! + . . \ Al (56) 

which is recognized as the series for r for a single slab where the exp —i — Di(x u —x L ) 
terms in A (41) have been expanded in a power series. 

3. List of Symbols 

E = electric wave vector 
H = magnetic wave vector 
H = static (earth's) magnetic field vector 
V= electron velocity vector 

jV= electron density, number of mediums, subscript for final medium 
m= electron mass 
v= electron collision frequency 



' When the limit of rf^O is taken Ai+x-^Ai and A 1+ i->A;. Also, as the matrices become diagonal the order of multiplication becomes immaterial. 
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e = dielectric constant for free space 
ju = permittivity for free space 

x\ 

V [= coordinate axes or variables of integration 

z 

a] 

„> = angles specifying direction of a vector in coordinate system 

Ph= direction of static field 

F=any wave field component 
71 [incident 

T\ = subscripts for I transmitted | quantities 
R\ [reflected J 

c = speed of light 

^ super- or subscript for ith layer, or for matrix quantities, or the imaginary unit 
m*— parameters appearing in quartic determinant — see (7) 

a\ 

dz,| = angular parameters for wave normal — see (8) 
a,\ 



parameters characterizing plasma — see (9) and (10) 



8) 

h\ 

h L \ 

hr) 

Dlf— variable in Booker quartic, root for 7th mode in ith layer 

/=wave frequency 
y = plasma frequency 
3h = gyro frequency 
V p = phase velocity 

co=wave angular frequency 

A= matrix relating field component of each mode to field components tangential to 

boundaries 
a^= components of A 

E\\ = component of E 7 lying in the plane containing the wave normal 
E±= component of E 7 lying perpendicular to plane containing the wave normal 

fti] 

; .r = ratios of electric field components for jth mode in ith layer — see (25) and (26) 
S}\ 

A t = matrix relating fields across ith layer — see (29) 

d t = thickness of iih layer 

M= matrix relating incident, transmitted, and reflected waves — see (32) 

E s , t, r,n, i, /,o= column matrices containing a component of E for each mode, see (20), (22), (27), 

and (34) 

N= matrix relating E s to £/ — see (37) 

II, ±R\\, jl = reflection coefficients 

T^= matrix product for ith layer — see (40) 

7' f ,"= matrix terms of series expansion of T t — see (44) 

1= identity matrix 

u, L= subscripts for upper and lower semi-infinite mediums respectively. 
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